using System;
using System.Collections;
using System.Collections.Generic;
using NUnit.Framework;
using dnAnalytics.Math;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Decomposition;
using dnAnalytics.LinearAlgebra.IO;
using dnAnalytics.Exceptions;

namespace dnAnalytics.UnitTests.LinearAlgebra.Decomposition
{
    public abstract class AbstractSvdTest
    {
        private static double ACCEPTABLE_ERROR = 1e-11;
        public abstract void GetDecompositionClass(String file, out Matrix matrix, out Svd svd);
        private Svd squareSvd;
        private Svd singularSvd;
        private Svd wideSvd;
        private Svd tallSvd;

        private Matrix squareMatrix;
        private Matrix singularMatrix;
        private Matrix wideMatrix;
        private Matrix tallMatrix;

        public abstract Vector GetVector(int size);
        public abstract Matrix GetMatrix(int order);
        public abstract Matrix GetMatrix(int rows, int columns);

        [TestFixtureSetUp]
        public void SetUp()
        {
            MatrixMarketReader mmr = new MatrixMarketReader();
            GetDecompositionClass("../../TestMatrices/random_real_general_array_100.matrix_market", out squareMatrix, out squareSvd);
            GetDecompositionClass("../../TestMatrices/gear_integer_general_coordinate_100.matrix_market", out singularMatrix, out singularSvd);
            GetDecompositionClass("../../TestMatrices/random_real_general_array_10_20.matrix_market", out wideMatrix, out wideSvd);
            GetDecompositionClass("../../TestMatrices/random_real_general_array_20_10.matrix_market", out tallMatrix, out tallSvd);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SvdConstructorNull()
        {
            Matrix matrix = null;
            Svd lu = new Svd(matrix, false);
        }

        [Test]
        public void Decomposition()
        {
            Matrix u = squareSvd.U();
            Matrix w = squareSvd.W();
            Matrix vt = squareSvd.VT();
            Matrix result = u * w * vt;

            double error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - squareMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - squareMatrix[i, j]) / squareMatrix[i, j]);
                    }
                    Assert.IsTrue(error < ACCEPTABLE_ERROR);
                }
            }

            u = tallSvd.U();
            w = tallSvd.W();
            vt = tallSvd.VT();
            result = u * w * vt;

            error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(tallMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - tallMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - tallMatrix[i, j]) / tallMatrix[i, j]);
                    }
                    Assert.IsTrue(error < ACCEPTABLE_ERROR);
                }
            }

            u = wideSvd.U();
            w = wideSvd.W();
            vt = wideSvd.VT();
            result = u * w * vt;

            error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - wideMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - wideMatrix[i, j]) / wideMatrix[i, j]);
                    }
                    Assert.IsTrue(error < ACCEPTABLE_ERROR);
                }
            }

            u = singularSvd.U();
            w = singularSvd.W();
            vt = singularSvd.VT();
            result = u * w * vt;

            error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - singularMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - singularMatrix[i, j]) / singularMatrix[i, j]);
                    }
                    Assert.IsTrue(error < ACCEPTABLE_ERROR);
                }
            }
        }

        [Test]
        public void DecompositionResults()
        {
            Matrix u = GetMatrix(squareMatrix.Rows);
            squareSvd.U(u);
            Matrix w = GetMatrix(squareMatrix.Rows, squareMatrix.Columns);
            squareSvd.W(w);
            Matrix vt = GetMatrix(squareMatrix.Columns);
            squareSvd.VT(vt);

            Matrix result = u * w * vt;

            double error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - squareMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - squareMatrix[i, j]) / squareMatrix[i, j]);
                    }

                    Assert.IsTrue(error < ACCEPTABLE_ERROR);
                }
            }

            u = GetMatrix(tallMatrix.Rows);
            tallSvd.U(u);
            w = GetMatrix(tallMatrix.Rows, tallMatrix.Columns);
            tallSvd.W(w);
            vt = GetMatrix(tallMatrix.Columns);
            tallSvd.VT(vt);

            result = u * w * vt;

            error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(tallMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - tallMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - tallMatrix[i, j]) / tallMatrix[i, j]);
                    }
                    Assert.IsTrue(error < ACCEPTABLE_ERROR);
                }
            }

            u = GetMatrix(wideMatrix.Rows);
            wideSvd.U(u);
            w = GetMatrix(wideMatrix.Rows, wideMatrix.Columns);
            wideSvd.W(w);

            vt = GetMatrix(wideMatrix.Columns);
            wideSvd.VT(vt);

            result = u * w * vt;

            error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - wideMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - wideMatrix[i, j]) / wideMatrix[i, j]);
                    }
                    Assert.IsTrue(error < ACCEPTABLE_ERROR);
                }
            }

            u = GetMatrix(singularMatrix.Rows);
            singularSvd.U(u);
            w = GetMatrix(singularMatrix.Rows, singularMatrix.Columns);
            singularSvd.W(w);
            vt = GetMatrix(singularMatrix.Columns);
            singularSvd.VT(vt);

            result = u * w * vt;

            error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - singularMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - singularMatrix[i, j]) / singularMatrix[i, j]);
                    }
                    Assert.IsTrue(error < ACCEPTABLE_ERROR);
                }
            }
        }

        [Test]
        public void DontComputeVectors()
        {
            Svd svd = new Svd(squareMatrix, false);
            Matrix w = svd.W();
            Assert.IsNull(svd.U());
            Assert.IsNull(svd.VT());

            Matrix result = GetMatrix(squareMatrix.Rows);
            svd.U(result);
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    Assert.AreEqual(0, result.ValueAt(i, j));
                }
            }

            svd.VT(result);
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    Assert.AreEqual(0, result.ValueAt(i, j));
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void UResultNull()
        {
            Matrix result = null;
            squareSvd.U(result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void VResultNull()
        {
            Matrix result = null;
            squareSvd.VT(result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void WResultNull()
        {
            Matrix result = null;
            squareSvd.W(result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SResultNull()
        {
            Vector result = null;
            squareSvd.S(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void UResultNotConformable1()
        {
            Matrix result = GetMatrix(squareMatrix.Rows + 1, squareMatrix.Rows);
            squareSvd.U(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void UResultNotConformable2()
        {
            Matrix result = GetMatrix(squareMatrix.Rows, squareMatrix.Rows + 1);
            squareSvd.U(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void VResultNotConformable1()
        {
            Matrix result = GetMatrix(squareMatrix.Columns + 1, squareMatrix.Columns);
            squareSvd.VT(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void VResultNotConformable2()
        {
            Matrix result = GetMatrix(squareMatrix.Columns, squareMatrix.Columns + 1);
            squareSvd.VT(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void WResultNotConformable1()
        {
            Matrix result = GetMatrix(squareMatrix.Rows + 1, squareMatrix.Columns);
            squareSvd.W(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void WResultNotConformable2()
        {
            Matrix result = GetMatrix(squareMatrix.Rows, squareMatrix.Columns + 1);
            squareSvd.W(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SResultNotConformable1()
        {
            Vector result = GetVector(squareMatrix.Columns + 1);
            squareSvd.S(result);
        }

        [Test]
        public void Rank()
        {
            int rank = squareSvd.Rank();
            Assert.AreEqual(squareMatrix.Rows, rank);

            rank = tallSvd.Rank();
            Assert.AreEqual(tallMatrix.Columns, rank);

            rank = wideSvd.Rank();
            Assert.AreEqual(wideMatrix.Rows, rank);

            rank = singularSvd.Rank();
            Assert.AreEqual(squareMatrix.Rows - 1, rank);
        }

        [Test]
        public void Norm()
        {
            double norm = squareSvd.Norm2();
            Assert.AreEqual(1011.4190766097751, norm, ACCEPTABLE_ERROR);

            norm = tallSvd.Norm2();
            Assert.AreEqual(445.27398550522514, norm, ACCEPTABLE_ERROR);

            norm = wideSvd.Norm2();
            Assert.AreEqual(385.51714390506197, norm, ACCEPTABLE_ERROR);

            norm = singularSvd.Norm2();
            Assert.AreEqual(2, norm, ACCEPTABLE_ERROR);
        }

        [Test]
        public void ConditionNumber()
        {
            double cn = squareSvd.ConditionNumber();
            Assert.AreEqual(207.5103850393096, cn, 1e-4);

            cn = tallSvd.ConditionNumber();
            Assert.AreEqual(5.01902319723716, cn, ACCEPTABLE_ERROR);

            cn = wideSvd.ConditionNumber();
            Assert.AreEqual(3.65923703921971, cn, ACCEPTABLE_ERROR);

            cn = singularSvd.ConditionNumber();
            Assert.Greater(cn, 1e16);
        }

        [Test]
        public void SolveSquareMatrix()
        {
            Matrix identity = MatrixBuilder.CreateIdentityMatrix(squareMatrix.Rows);
            Matrix result = squareSvd.Solve(identity);
            result = squareMatrix * result;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(1, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                }
            }
        }

        [Test]
        public void SolveTallMatrix()
        {
            Matrix matrix = GetMatrix(tallMatrix.Rows, 5);
            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    matrix[i, j] = (i + 1.0) / (j + 2.0);
                }
            }
            Matrix result = tallSvd.Solve(matrix);
            Matrix test = (tallMatrix.Transpose() * tallMatrix).Inverse() * tallMatrix.Transpose() * matrix;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    Assert.AreEqual(test[i, j], result[i, j], Constants.ACCEPTABLE_ERROR);
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNull()
        {
            Matrix matrix = null;
            squareSvd.Solve(matrix);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformable()
        {
            Matrix matrix = GetMatrix(squareMatrix.Rows + 1);
            squareSvd.Solve(matrix);
        }

        [ExpectedException(typeof(InvalidOperationException))]
        public void SolveMatrixInvalidOperation()
        {
            Svd svd = new Svd(squareMatrix, false);
            Matrix matrix = GetMatrix(squareMatrix.Rows);
            svd.Solve(matrix);
        }

        [Test]
        public void SolveSquareMatrixResult()
        {
            Matrix identity = MatrixBuilder.CreateIdentityMatrix(squareMatrix.Rows);
            Matrix result = GetMatrix(squareMatrix.Rows);
            squareSvd.Solve(identity, result);
            result = squareMatrix * result;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(1, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                }
            }
        }

        [Test]
        public void SolveTallMatrixResult()
        {
            Matrix matrix = GetMatrix(tallMatrix.Rows, 5);
            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    matrix[i, j] = (i + 1.0) / (j + 2.0);
                }
            }
            Matrix result = GetMatrix(tallMatrix.Columns, 5);
            tallSvd.Solve(matrix, result);
            Matrix test = (tallMatrix.Transpose() * tallMatrix).Inverse() * tallMatrix.Transpose() * matrix;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    Assert.AreEqual(test[i, j], result[i, j], Constants.ACCEPTABLE_ERROR);
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNullResult1()
        {
            Matrix matrix = null;
            Matrix result = GetMatrix(squareMatrix.Rows);
            squareSvd.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNullResult2()
        {
            Matrix matrix = GetMatrix(squareMatrix.Rows);
            Matrix result = null;
            squareSvd.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformableResult1()
        {
            Matrix matrix = GetMatrix(squareMatrix.Rows + 1);
            Matrix result = GetMatrix(squareMatrix.Rows);
            squareSvd.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformableResult2()
        {
            Matrix matrix = GetMatrix(squareMatrix.Rows);
            Matrix result = GetMatrix(squareMatrix.Rows + 1);
            squareSvd.Solve(matrix, result);
        }

        [ExpectedException(typeof(InvalidOperationException))]
        public void SolveMatrixInvalidOperationResult()
        {
            Svd svd = new Svd(squareMatrix, false);
            Matrix matrix = GetMatrix(squareMatrix.Rows);
            Matrix result = GetMatrix(squareMatrix.Rows);
            svd.Solve(matrix, result);
        }

        [Test]
        public void SolveSquareVector()
        {
            Vector vector = GetVector(squareMatrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (i + 1.0) / (i + 2.0);
            }
            Vector result = squareSvd.Solve(vector);
            result = squareMatrix * result;

            for (int i = 0; i < result.Count; i++)
            {
                Assert.AreEqual((i + 1.0) / (i + 2.0), result[i], Constants.ACCEPTABLE_ERROR);
            }
        }

        [Test]
        public void SolveTallVector()
        {
            Vector vector = GetVector(tallMatrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (i + 1.0) / (i + 2.0);
            }
            Vector result = tallSvd.Solve(vector);
            Vector test = (tallMatrix.Transpose() * tallMatrix).Inverse() * tallMatrix.Transpose() * vector;

            for (int i = 0; i < result.Count; i++)
            {
                Assert.AreEqual(test[i], result[i], Constants.ACCEPTABLE_ERROR);
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNull()
        {
            Vector vector = null;
            squareSvd.Solve(vector);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformable()
        {
            Vector vector = GetVector(squareMatrix.Rows + 1);
            squareSvd.Solve(vector);
        }

        [ExpectedException(typeof(InvalidOperationException))]
        public void SolveVectorInvalidOperation()
        {
            Svd svd = new Svd(squareMatrix, false);
            Vector vector = GetVector(squareMatrix.Rows);
            svd.Solve(vector);
        }

        [Test]
        public void SolveSquareVectorResult()
        {
            Vector vector = GetVector(squareMatrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (i + 1.0) / (i + 2.0);
            }
            Vector result = GetVector(squareMatrix.Rows);
            squareSvd.Solve(vector, result);
            result = squareMatrix * result;

            for (int i = 0; i < result.Count; i++)
            {
                Assert.AreEqual((i + 1.0) / (i + 2.0), result[i], Constants.ACCEPTABLE_ERROR);
            }
        }

        [Test]
        public void SolveTallVectorResult()
        {
            Vector vector = GetVector(tallMatrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (i + 1.0) / (i + 2.0);
            }
            Vector result = GetVector(tallMatrix.Columns);
            tallSvd.Solve(vector, result);
            Vector test = (tallMatrix.Transpose() * tallMatrix).Inverse() * tallMatrix.Transpose() * vector;

            for (int i = 0; i < result.Count; i++)
            {
                Assert.AreEqual(test[i], result[i], Constants.ACCEPTABLE_ERROR);
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNullResult1()
        {
            Vector vector = null;
            Vector result = GetVector(squareMatrix.Rows);
            squareSvd.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNullResult2()
        {
            Vector vector = GetVector(squareMatrix.Rows);
            Vector result = null;
            squareSvd.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformableResult1()
        {
            Vector vector = GetVector(squareMatrix.Rows + 1);
            Vector result = GetVector(squareMatrix.Rows);
            squareSvd.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformableResult2()
        {
            Vector vector = GetVector(squareMatrix.Rows);
            Vector result = GetVector(squareMatrix.Rows + 1);
            squareSvd.Solve(vector, result);
        }

        [ExpectedException(typeof(InvalidOperationException))]
        public void SolveVectorInvalidOperationResult()
        {
            Svd svd = new Svd(squareMatrix, false);
            Vector vector = GetVector(squareMatrix.Rows);
            Vector result = GetVector(squareMatrix.Rows);
            svd.Solve(vector, result);
        }

    }
}
